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ABSTRACT 

Using a one level, barotropic ocean model, driven by sur- 
face winds, a finite difference form of the vorticity equation 
was integrated over 210 days of simulated time. The solutions 
using constant coefficients of lateral eddy viscosity were 
compared with those using variable coefficients derived from 
enstrophy cascade and energy cascade. Using a constant eddy 
viscosity coefficient of rather low magnitude produces a large 
amplitude computational oscillation which fills the entire 
basin. An order of magnitude larger coefficient produces a 
marginally satisfactory solution, where the western boundary 
current was rather well represented, but a moderate computa- 
tional oscillation was still evident. By increasing the co- 
efficient yet another order of magnitude, the computational 
oscillation is negligible, but the solution in the ocean in- 
terior is unrealistically damped. An accurate physical and 
numerical depiction of both the ocean interior and western 
boundary with no computational oscillation was achieved by 
using either of the two forms of non linear eddy viscosity. 
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I. 



INTRODUCTION 



This study is part of a broader effort to develop the 
capability of making large scale oceanographic predictions 
on a dynamical basis. The large-scale and meso-scale thermal 
structure of the ocean is very dynamic and is related to most 
marine and atmospheric processes. Some of the most obvious 
relationships to sea surface temperature (SST) are ocean 
fronts, circulation, currents, and sea state; and atmospheric 
temperature, circulation, and wind velocity. The effect that 
these environmental phenomena can have on naval and maritime 
operations is well known. 

The significance of the long range effect of SST anomalies 
on weather patterns has received increasing scientific aware- 
ness. On 20 August 1974, J. Namias, at a NORPAX Conference 
at the U.S. Naval Postgraduate School, reported on the results 
of an empirical ocean/atmosphere prediction model. With a 
knowledge of SST anomalies in the spring of 1974, Namias 
derived fields of atmospheric temperature anomalies and cir- 
culation for the following summer. Namias predicted in May 
that a comparatively severe drought would occur over the 
midwestern United States during the summer of 1974. Later 
events verified his forecast. 

One of the areas of difficulty in adequately representing 
the circulation and anomalies in a finite-grid ocean circula- 
tion model, and subsequently satisfying and integrating the 
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equations of motion, is the representation of internal 
frictional effects within the ocean fluid. 

Using linear theory, Takano [1974] showed that if the 
coefficient of eddy viscosity is too small, a computational 
oscillation in space results from not resolving the western 
boundary current. This computational oscillation fills the 
entire basin and contaminates the solution. On the other 
hand, if the coefficient is increased, the open ocean solu- 
tion is too viscous. To handle this problem, Takano showed 
that the use of upstream differencing in the beta term of 
the vorticity equation allows a smaller coefficient than per 
mitted by linear theory. However, this method unfortunately 
produces excessive damping of time dependent motion. In 
addition, this mathematical scheme is not representative of 
any physical motion in the ocean. A better approach, from 
a physical point of view, is the use of non linear eddy vis- 
cosity. The friction force which arises using non linear 
eddy viscosity is extremely sensitive to the scale of the 
motion, and therefore will be relatively small in the oceans 
interior where the scales are comparatively large, and will 
be relatively large in the western boundary region where the 
scales are comparatively small. The comparatively large 
dissipation in the western boundary region will keep the 
current broad enough to be resolved by the grid and thereby 
prevent the formation of any computational oscillation. 

In 1968, Leith examined two dimensional turbulence 
advection and derived non linear coefficients of eddy 
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viscosity from the cascade of enstrophy from large scale to 
small scale motion. In this case, the coefficient of eddy 
viscosity is proportional to the magnitude of the horizontal 
gradient of vorticity. 

Another satisfactory non linear procedure was introduced 
by Smagorinsky [1963] in which an estimate of the energy cas- 
cade rate is made from the fluid deformation itself. In this 
case, the coefficient of eddy viscosity is proportional to 
the absolute value of the deformation. 

The purpose of this thesis is to introduce non linear 
coefficients of lateral eddy viscosity that are not only 
reasonably representative of actual physical conditions, but 
these same coefficients must retain their reasonable repre- 
sentation when used in finite grid numerical ocean circulation 
models. In addition, with these coefficients the numerical 
model must remain computationally stable in long term inte- 
gration. Integrations using both of the above forms of non 
linear eddy viscosity are examined and compared with numerical 
and analytic results using constant coefficients. 
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II. THE MATHEMATICAL STATEMENT OF THE PROBLEM 



A. FORM OF THE VORTICITY EQUATION 

The continuity equation for a homogeneous fluid and the 
non linear equations of horizontal motion, cross differenti- 
ated to eliminate the pressure term, form the basis for the 
vorticity equation: 

dr 3 T y ^ T x 

dt + v 3 + fV- V = 3"x(F y + p^T) - 9y(F x + p-^H) 

or |^+V-V^+vB + fv-V = ^C F y + IX_J • *|:( F X + Tx ) 

9t PoH dy poll 

Cil-i) 

A -F 

where 3 = is taken constant. In (II-l), ( T x, T y) is the 
surface stress, H is the constant basin depth, and all other 
symbols have their usual meaning. Bottom stress was neglected. 
The friction forces are represented by 
F x = V- (Ay u) 

F y = V • (A v v) (II-2) 

where A is the coefficient of lateral eddy viscosity and (u,v) 
is the vertically averaged velocity, as discussed below. 

The model is designed with a ’’rigid surface,” in order 
to filter gravity waves from the system: 
w(o) = 0. 

It is also necessary that there be no vertical motion on the 
flat ocean floor: 
w(-H) = 0. 

Because w is zero at the top and bottom, the divergence of 
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the vertically averaged current, V, is zero. 
V • V = 0. 



Therefore V can be defined by a streamfunction , if; , 
u - - 

3y 

7 = life 

3X. 

Integrating (II-l) from top to bottom and using y V = 0, 
the final form of the vertically integrated vorticity equation 
was written 



at * 6 H- ■ hC' v ' 7v * p y + ^ 



- |yC- V- VU ♦ F x - T ^) 

Equation (II-3) also includes the assumption that 



(I I ~ 3) 



V • Vv = V • v , e tc . 

The next section will discuss the formulation of the 
friction terms, F x and Fy, which depend upon the non linear 
eddy viscosity coefficient, A. 



B.. ENSTROPHY CASCADE FORM OF NON LINEAR EDDY VISCOSITY 



The first method of generating non linear coefficients 
of eddy viscosity is through the theory of two dimensional 
turbulence, in which enstrophy and kinetic energy are con- 
served by the advective terms. Using these principles, 
Kraichman [1967] and Leith [1968] derived the relation 

E (k) ec ri 2 / 3 k ” 3 (II -4) 

where q is the assumed constant cascading rate of mean squared 
vorticity and E(k) is the kinetic energy in wave number k. 

The eddy viscosity which causes the dissipation is also assumed 
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a function of n and k, and by dimensional analysis 

A = a n V 3 k" 2 * ( 1 1 - 5) 

where a is a constant. 

One estimate of n made locally as a function of sur- 
rounding data was made by Leith 

n = A( V O ’ ( V O = A | V£ | 2 , (II-6) 

where 5 is the local vertical component of vorticity. 
Substituting (II - 6 ) into ( 1 1 - 5 ) and solving for A leads to 

A = ( « V*/k) 3 I VC | . 

Assuming that the wave number k* = 2u/d, where d is the grid 
size, lies in the inertial subrange, the final form of non 
linear eddy viscosity for enstrophy cascade becomes 

A = (a 1 / 1 3 1 V? | . CII-7) 

Denoting the minimum viscosity coefficient by A Q , the 
viscosity in the model was written, 

A = A 0 + y | V£ |d 3 (I I - 8) 

A is a maximum when | V£ | is maximum. Defining A max = mA 0 , 
then (II-8) becomes 

•A-max = ^o = A 0 + 1 1 Imax ^ 

Solving for y 

_ Ao(m-l) 

Y “I 71 [ cF 

1 ^ max a (II - 10) 

The quantity | V c | max was estimated by 

I I max = 6.67x10 cm' 1 sec' 1 (11-11) 

which for d = 300 km leads to a value for y of 
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Y 



ApO-1) 

1.8 x 10 9 



QI-12) 



Substituting this into (II -8) gives the final equation for 
the enstrophy cascade form of non linear eddy viscosity for 
this study: 

A = A 0 (1 + = A 0 (l + Cm - l)i v d 0 (11-13) 

I v ? I max 6.67x10 

In (11-13), H' inax was taken from linear theory [Munk, 1950], 
and m was considered an adjustable parameter. 

With this non linear form of eddy viscosity, the vertically 
averaged friction force terms in the x and y directions become 
from (I I -2) 



= -i-( A^ii) + 


JL( aL“) 


(11-14) 


9x 3x 


3y ay 




- -,-2-C A f) * 


J-ca|1) . 


(11-15) 


ax a x 


ay a y 





The vertical component of the curl of the friction force 
becomes 

CURL z F = ilx _ a _lx . (11-16) 

ax a y 

C. KINETIC ENERGY CASCADE FORM OF NON LINEAR EDDY VISCOSITY 
From diminsional arguments, Leith [1968] showed that in 
three dimensional turbulence, if the energy cascade rate from 
large scale to small scale is proportional only to kinetic 
energy dissipation (e) and wave number (k) , then 

E (k) = a x e 2/3 k" 5 / 3 . (11-17) 

If, in addition to this equation, it is assumed that eddy 

viscosity, which produces the dissipation, is also a function 

of e and k only, then by dimensional arguments 

15 



A = a 2 e V 3 k ” V 3 (11-18) 

where k lies in an inertial subrange. Assuming the dissipa- 
tion rate is constant, then k ~ 4 / 3 % d V 3 . This quasi-linear 
eddy viscosity is dependent only on grid size, and hence 



latitude if the grid size is 
realistically, dissipation is 
of motion or grid size. Smag 
value of dissipation: 
c = A | D | 2 

where | D| is the magnitude of 
this case (11-18) becomes 

A = a 2 3 / 2 | D | d 2 

In the deformation tensor |D| 

deformation is Do = £v. 

5 Jx dy 

Dt = 3u - |v . 

3 x Jy 

The Smagorinsky equation 
( 1 1 - 8 ) for this study 
A = Ao + Y I D | d 2 . 



atitude dependent. But more 
not constant nor independent 
rinsky [1963] assumes a local 

(11-19) 

the deformation tensor. In 

( 11 - 20 ) 

= / D s z + D t i! , the shearing 
and stretching deformation is 

akes the form similar to 

(II-21) 



Defining A max = mA 0 , equation (11-21) becomes 

mA 0 - A 0 + Y | D | max d 

Solving for y 

= Ao(m-l) 

I max 



( 11 - 22 ) 



(11-23) 



Using (11-23) in 11-21), the final energy cascade form of 
the non linear eddy viscosity becomes 
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( 11 - 24 ) 



(m-1) 


D 




D 



] 

max 

The form of the friction force principally used in the 
model for the kinetic energy cascade case was of the form 
of (11-14) and (11-15). Experiments following the friction 
force of Smagorinsky [1963] were also used for comparison: 



F x = -2- (A D t ) + (A D s ) 

3x 3y 

F y = _2_ (A D s ) - -2- (A D t ) 



(11-25) 

(11-26) 



3x " sy 

In both cases the curl of the friction force term was (11-16) 
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III. THE MODEL AND FINITE DIFFERENCE EQUATIONS 



A. PHYSICAL CHARACTERISTICS OF THE MODEL 

The model consisted of a one-level barotropic ocean in a 
square basin of length, L = 9600 km; of breadth, B = 9600 km; 
and a flat bottom of constant depth, H = 2 km. A portion of 
the staggered grid is shown in Figure 1. The total grid 
points (excluding the boundary buffer) are 33 x 33, and the 
distance between adjacent grid points is 

X = Y = d = 300 km. (I II - 1) 

The value chosen for g corresponds to a grid centered at 
32.5° N. 

The staggered grid consists primarily of: 1) the inter- 
sections, or principal grid points ( • ), where are defined 
the streamfunction ( T ) , vorticity ( r; ) , deformation ( D ) , 
and the coefficient of eddy viscosity (Ap) generated from 
deformation; and 2) the grid centers (o), where are defined 
velocity (u,v) , the gradient of vorticity the friction 

force (F x , Fy) , and the coefficient of eddy viscosity (A^ ) 
generated from (Vt;). Auxiliary variables were defined at 
cross ( x ) points . 

The model is not strictly designed to simulate any partic- 
ular ocean, but rather to be representative of the fundamental 
physical characteristics of mass transport and western boundary 
current of an ocean in the northern hemisphere as affected by 
the viscous action of the ocean's large scale circulation. 
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i 



l 





2,1 



J. 




Figure 1: The Staggered Grid with Corner Boundary. The solid line represents 

the physical boundary of the basin* 

(• ) Principal grid point ( 0 , £ , D, A( |]}| ) ) 

( o) Grid center (u, v, , A( V ), F , F ) 

x y 

(x) Grid average points (UAV, VAV, UAVE, VAVE, A(|D|) ave , , F ygve ) 
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The energy for this barotropic wind driven circulation 
model is represented by 



'o 



= - F cos(^), T y = 0, (HI - 2) 

a pattern of westerly winds in the center half and easterly 
winds in the northern and southern quarters of the ocean 
(See Figure 2) . This leads to a stress curl term 



. 9 



t_ r I2L_ > _-2nf . , 

°y''PoH-)j HB sin ( 



2ny 



A 



x 



(III -3) 



B ' 



which provides the actual forcing in the vorticity equation. 
F is the amplitude of the zonal component of the wind, and 
B is the north-south extent of the domain. 

The velocity (u,v) written in numerical form is 



Uij = jj-K n-i.i-i + Ti t j-i ) - ( Ti-ij + nj )3 

2 2 

Vi.i = r c^i.3 + > j _ 1 C^i" 1 » j + yi-l.J-1)] 

J 2 2 

The boundary conditions for (u,v) used along the coastline 

were zero normal flow and zero slip. Zero normal flow was 

attained by requiring the streamfunction ( Y ) on the left side 

of equation ( 1 1 - 3) to be equal zero on the ocean perimeter. 

To implement the condition of zero normal flow and zero slip 

in the terms on the right hand side of (11-33 > the velocity 

is defined as zero on the coastline by defining the zonal 

boundaries 

u i , 1 = -u i>2 u i , 34 = -ui,33 

vi f l = -ui } 2 Vi } 34 = -vi,33, (HI-5) 

for the meridional boundaries 
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Figure 2: Wind Stress Pattern 
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u l,j = _u 2 , j 



u 34,j - ~ u 33 , j 



and for the corners 

(u > v) 1 ^ 1 = (u,v ) 2>2 

Cu, v) i ^ 24 = (u,v ) 2;33 



v 34,j - ' v 33,j, 
Cu,v ) 34;1 = Cu,v) 33 ;2 

( u » v ) 34 f 34 - Cu,v) 33> 33 



(1 1 1 - 6 ) 



(III-7) 



where the original 32 x 32 (u,v) grid field becomes a 
34 x 34 grid to include the boundary conditions. 



B. FINITE DIFFERENCE FORM OF THE VORTICITY EQUATION 
The finite difference form of (II- 3) was written 

JL y2 u/ . . +6 ( y i + i>3 ~ _ 

3t i,l P 2Ax 

1 (V* Vv)^ +1 -j + i + (V*Vv)i + 1 

’ TIC 2 ) 

- (1 V * vv )j, i + 1 + CV-Vv)i t i }] 

i rr CV-vu ) ijj+1 + cv-vu) i+1)j+1 ^ 

d U 2 J 

CV-VujjJ + (V-Vu ) i+1>1 



+ i [C p yi + i , ULi . ( F yj,i + i + Fyi »j )] 

- -*• r r px i.i + l + px i+l.i + l -) . r px i.i ^ Px i + l,j 'v i 
<1 1 K 2 2 



■ ttt sin y P d 

In (III- 8 ) V 2 Y . . is given by 



( 1 1 1 - 8 ) 



V 2 y . 



1,1 



p + ♦ * 
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i > j + 1 * 



4 Vii ] - 



The advection terms ( V-VU, etc.) were finite differenced 
using the same method as Haney [1974], and the friction 
terms were expressed differently for each of the two forms 
of eddy viscosity as discussed below. Readers are referred 
to the documented computer program in Appendix A for further 
details . 

In order to define F, it is first necessary to determine 
three values of minimum A (AMIN) , such that the maximum in 
the analytic s treamfunction (T) is at x = d/2, d, 2d, 
respectively. Since the grid cannot resolve the western 
boundary in the first case (where the western boundary width 
is d) , a large computational mode was expected with AMIN(l). 
The second case was expected to be marginally stable with 
AMIN(2); and finally, with AMIN(3) it was expected that the 
western boundary would be clearly resolved, but the value of 
AMIN(3) would be so high as to unrealistically damp the 
interior ocean solution. The respective difference equations 
for the three AMINs were: 



AMIN (1) = 


(3 


x {l 
2n 


x d/2) 3 




AMIN (2) = 


(3 


* 7 1 
2n 


x d) 3 




AMIN (3) = 


(3 


X {1 


x 2d) 3 


(III-9) 




2n 





These three values of AMIN were used for the three fundamental 
experiments where A was taken as constant, and also for the 
minimum A in the experiments with non linear coefficients of 
eddy viscosity using (11-13) and (11-24) with A 0 = AMIN. 
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In determining the grid size and location for the non 
linear coefficients of eddy viscosity, the grid location of 
the generating parameters ( | V£ | or |D|) had to be taken into 
consideration. This led to a 34 x 34 grid for A(|V?|) in 
the case of enstrophy cascade, and a 33 x 33 grid for A(|D|) 
in the kinetic energy cascade form. (See Figure 1) . 

1 . Enstrophy Cascade Case 

In the enstrophy cascade case, in order to generate 
coefficients of non linear eddy viscosity, it is necessary 
to generate a relative vorticity field and to determine the 
gradient of vorticity at each time step. Several possible 
techniques exist in developing a vorticity field, where in 
all cases vorticity would be defined on the 33 x 33 principal 
grid points (•)• For this model C was defined in terms of 
(u,v) by 



C = 



3v 

3x 



3u 

3y 

VAV 



'I >J 



i. +AJ - 



VAV i , j . UAV i,H - UAV i,,j 



Ci i i - io) 



i = 1...33, j = 1. .33 
where UAV and VAV are defined as 



UAV • . = Ul >J + Ui+1 >J ; VAVi A 
1 > J 2 1 > J 



Vi, 3 + 



i = 1...33, j = 1...34; i = 1...34, j 






(III-ll) 

= 1 , «33 



in order to get an equivalent (u,v) on the principal grid 
points. If (u,v) are written in terms of V, then (III-10) 
reduces to 



= 2~p ( A i + 1 J + 1 + f i+l,j-l 



Yi-lJ + 1 - 4 n,j) 
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+ "i-U-1 



(1 1 1 - 12 ) 



This method, along with three other schemes of Miyakoda [1962] 
for calculating vorticity, were utilized and compared. 

The gradient of vorticity (|vc|) was developed using 
centered differences on the 32 x 32 grid centers ( ° ) : 



1 vc I i , j { [ i r c 1 , 1 * , i -1 ■ * s i-l,i ]» 






+ ‘’l' 



% ' S A.J_ + . I hlA il ) ] 2 > ^ 

2 2 



i = 2. . .33, j = 2, . .33 



(1 1 1 - 1 3 ) 



Finally, (11-13) in finite difference form 
A(|Vc|)ij = AMIN x (1 + m x | I i s j 

I ^ I max (1 1 1 - 14) 

i = 2. ..33, j = 2. ..33 

where Ai } j is defined on the 32 x 32 grid centers. The value 
of | V? | m of 6 x 10 1 4 was estimated from linear theory 
[Munk, 1950], and m was considered an adjustable parameter. 

The finite difference form of the friction force for 
enstrophy cascade, on a 32 x 32 grid, was written from (11-14) 
and (11-15) as 

Fx i ■ = ~ [( A i ? j + Ai+ l ? l)( u i + ;L,j - u i,j) - 
> j d 2 2 

- ( A i-1.1 + A i > 1 ) ( U i » J ‘ U i' 1 > j ) ] 

2 

+ !-> [( Ai »-i + A i,j + l ){ u i,i + l ■ u i,j) 

- r A i, j-1 + A i,j ) ( u i» 3 " U iJ-l)] 

2 



1 = 2. . .33, j = 2. . .33 
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= k ^-; Ai - +i, j )(vui,j - vi, j 



j) 



- ( Ai : . 1 > j . I A i . J )C v i,j - v i _ 1 » j ) 1 



+ \ 2 ( Vi > J +1 ■ v nj) 

- C A i > 3 - 1 + A i,i )(vi,j - v i , j - 1) ] (I II - 15) 

2 

i = 2. ..33, j = 2. ..33 

It is clear that (III-15) requires a 34 x 34 field 
of A(|V£|) in order to calculate the friction forces. Two 
methods were utilized to generate a value of |V^| and thereby 
a value of eddy viscosity coefficients across the boundary. 

The first case was linear extrapolation in all four direc- 
tions. For the western boundary the finite difference form 
of A( | v? | ) was 

SLOPE = (A, ■ - A, H/d 

A, = SLOPE x d + A~ (1 1 1 - 16) 

j- y J L. y J 

j = 2 , ... 33 

and the values of A immediately outside the other boundaries 
were determined in a similar manner. This method is physically 
representative because it gives a boundary of increasing co- 
efficients in the direction of the western boundary. Along 
the other boundaries the gradient, and hence the boundary 
value of A(|Vs|) is flat. 

The second method utilized to generate a value for 
A(|V?|) across the boundary was to extend the existing interior 
"boundary" outwards in all four directions, so that across the 
western boundary, for example 
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A 



1> j 



l 2, j 



j = 2 33 



CHI -17) 



This method is consistent with the antisymmetric velocity 
profile which accompanies the zero flow boundary conditions. 
In both cases the 34 x 34 corners were not used for A(|Vt|). 

The curl of the friction force, defined on a 31 x 31 
interior grid, was 



of 



CURL F i} j = ( Fy i ,j + p yi,j-l - Fy i-l,j + _ jXj H , j ~ 1 )/d 

2 2 

. ( Fx i-l.i + Fx i.j ' Fx i-l,j-l + Fx i > j - 1 j /d 
2 2 

i = 3 33 j = 3 33 (III-18) 

2 . Kinetic Energy Cascade Case 

The first steps in determining the deformation field 
the fluid was 1) to calculate the shearing deformation 



(D s ) 

D , VAVj.l.j ' VAVj.j t UAVj.j.l - UAVj.j 

’J d d (I II - 19) 

i = 1. .33 j = 1 .... 33 

where UAV and VAV are defined in (III-ll), and 2) to calculate 
the stretching deformation (D^-) 



D t . . = UAVE i+1 j - UAVE i} j VAVEij+2 - VAVE i} j 
1,J d ’ d 



(III -20) 



i = 1. .33 j = 1 .... 33 
where UAVE and VAVE are defined as 

i = 1 , . . .34 
j = 1 , ... 33 



UAVE i ■ = U i,i + U i , j + 1 
X > J 



VAVEi j = V L , . 1 . + V i-l, .j 

2 



i = 1 .... 33 
j = 1 34 



(III -21) 



Then the square root of the squares of the values of 



27 



deformation along and normal to a streamline gave the value 
of deformation at each principal grid point (•): 



|D|i,j = [( Ds i,j) 2 + ( D t i j) 2 ] 1 / 2 (III-22) 

i = 1 .... 33 j = 1 .... 33 

Next, equation (11-25) in finite difference form was used to 
calculate values for A(|d|) based on kinetic energy cascade 

A C I D | ) - ^ = AMIN(1 - m x iEijji.) (III-23) 

1,3 | D max | 

i = 1 .... 33 j = 1 .... 33 

where AMIN and m are as described above, and the maximum 
deformation (D max ) was estimated to be 6.22 x 10” 7 . Since 
A(|D|) was a 33 x 33 field, the form of the friction force 
did not require additional boundary conditions for A(|D|). 

The numerical form for the friction forces based on 
A(|D|) written from (11-14) and (11-15) was 



x i > j ' ^ [( 



1 A-; -j _ -j Aj -j 

,J 1 — )( u i + l,j - u i,j) 



d‘ •• 2 

- ( A i ~ 1 > 3 - 1 + A i -1 > .1 ) ( u i > 3 - u i-l,j)] 

2 

+ a 2 f^ Ai - »i + ^ i ' 1 > .i) ( U iJ + 1 - U iJ) 

2 

- ( A i>j + A i-l>j-l ) r u i,j - u i,j-l)] 



F Yi , j 



= T2 + A i>j ) ( v i + l , j - v i,j 



' 3 ) 



- + A i-l,3 H v i>i ' Vi - 1 >3)] 



* n 2 U A i,j * A i-l,i)t v i.i*l - v i > j ) 
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